# Modelo de hiv em Banks et al, 2012

params <- c(
	aA  = 1e-12,
	aE  = 0.1,
	aT  = 0.008,
	bE1 = 0.328,
	bE2 = 0.001,
	bT1 = 1e-2,
	bT2 = 1e-6,
	cE  = 0.001,
	cT  = 7,
	cV  = 13,    # natural clearance rate of Vi
	dmE = 0.25,  # maximum death rate of E1 due to infected Th
	drE = 0.01,  # rate at which E1 eliminates T1*
	dE1 = 0.1,
	dE2 = 0.005,
	dT1 = 0.02,
	dT2 = 0.005,
	dV  = 0.7,
	gE  = 0.01,
	gT  = 0.005,
	nE  = 3,
	nT  = 2,
	nV  = 100,
	kV  = 0.1,
	kb1 = 0.1,
	kb2 = 100,
	kd  = 0.5,
	kg  = 10,
	ks  = 1e2,
	sE  = 0.001,
	nu  = 10  # volume scale
)

a <- c(
	"dT1*x1",
	"drE*(x6/nu)*x2",
	"dT2*x3",
	"dT2*x4",
	"cV*x5",
	"(dE1+dmE*(x2/(x2+kd)))*x6",
	"dE2*x7",
	#####################
	"bT1*(x5/nu)*x1",
	"bT2*(x5/nu)*x3",
	#####################
	"gT*x1",
	"gT*x2",
	"gE*((x1+x2)/(x1+x2+kg))*x6",
	#####################
	"(aT*(x5/(x5+kV))+aA)*x3",
	"(aT*(x5/(x5+kV))+aA)*x4",
	"aE*(x5/(x5+kV))*x7",
	#####################
	"dV*x2",
	"(cT*(ks/(x5+ks)))*nu",
	"nu*cE+bE1*(x2/(x2+kb1))*x6",
	"bE2*(kb2/(x7+kb2))*x7"
)

V <- t(matrix(c(
	-1, 0, 0, 0, 0, 0, 0,
	 0,-1, 0, 0, 0, 0, 0,
	 0, 0,-1, 0, 0, 0, 0,
	 0, 0, 0,-1, 0, 0, 0,
	 0, 0, 0, 0,-1, 0, 0,
	 0, 0, 0, 0, 0,-1, 0,
	 0, 0, 0, 0, 0, 0,-1,
	#####################
	-1,+1, 0, 0,-1, 0, 0,
	 0, 0,-1,+1,-1, 0, 0,
	#####################
	-1, 0,+1, 0, 0, 0, 0,
	 0,-1, 0,+1, 0, 0, 0,
	 0, 0, 0, 0, 0,-1,+1,
	#####################
	params['nT']*+1,0,-1,0,0,0,0,
	0,params['nT']*+1,0,-1,0,0,0,
	0,0,0,0,0,params['nE']*+1,-1,
	#####################
	0,-1,0,0,params['nV']*+1,0,0,
	0, 0,+1, 0, 0, 0, 0,
	0, 0, 0, 0, 0,+1, 0,
	0, 0, 0, 0, 0, 0,+1
), ncol=7, byrow=TRUE))

X0 <- params['nu'] * c(x1=5, x2=1, x3=1400, x4=1, x5=10, x6=5, x7=1)

# algorithm

library(GillespieSSA)
X <- ssa(X0, a, V, params, 100, 'BTL')

# plot

library(gridExtra)
library(ggplot2)
library(reshape2)

theme_set(theme_gray(base_size=10) %+replace% theme(
	axis.text        = element_text(color='black'),
	panel.grid.minor = element_line(color='gray80'),
	panel.grid.major = element_line(color='gray50'),
	panel.background = element_rect(fill='white'),
	legend.text.align=0,
	legend.key = element_blank()
))

df <- as.data.frame(X$data)
colnames(df)[1] <- 'time'
rownames(df) <- 1:nrow(df)
df <- melt(df, id='time', variable.name='type', value.name='counts')

labels <- c(x1=expression(CD4[act]), x2=expression(CD4[act]^'\u2217'), x3=expression(CD4[rest]), x4=expression(CD4[rest]^'\u2217'), x5=expression(V[I]), x6=expression(CD8[eff]), x7=expression(CD8[mem]))

p <- ggplot(df, aes(x=time, y=counts, color=type, size=type, linetype=type)) + geom_line()
p <- p + scale_y_log10()
p <- p + scale_colour_manual(name="Type", values=c(x1='blue',x2='blue',x3='blue',x4='blue', x5='red', x6='green',x7='green'), labels=labels)
p <- p + scale_size_manual(name="Type", values=c(x1=1,x2=1,x3=.5,x4=.5,x5=.5,x6=1,x7=.5), labels=labels)
p <- p + scale_linetype_manual(name="Type", values=c(x1=1,x2=2,x3=1,x4=2,x5=1,x6=1,x7=1), labels=labels)
print(p)
ggsave('plot.pdf', p)



## ODE

names(X0) <- c('T1','T1_','T2','T2_','VI','E1','E2')

fn_ode <- function(t, X, params) {
	with(as.list(c(X, params)), {
		# I use a double 'dd' prefix to avoid name collision
		ddT1  <- - dT1*T1 - bT1*VI*T1 - gT*T1 + nT*((aT*VI)/(VI+kV)+aA)*T2
		ddT1_ <- bT1*VI*T1 - dV*T1_ - drE*E1*T1_ - gT*T1_ + nT*((aT*VI)/(VI+kV)+aA)*T2_
		ddT2  <- cT*(ks/(VI+ks)) + gT*T1 - dT2*T2 - bT2*VI*T2 - ((aT*VI)/(VI+kV)+aA)*T2
		ddT2_ <- gT*T1_ + bT2*VI*T2 - dT2*T2_ - ((aT*VI)/(VI+kV)+aA)*T2_
		ddVI  <- nV*dV*T1_ - cV*VI - (bT1*T1 + bT2*T2)*VI
		ddE1  <- cE + ((bE1*T1_)/(T1_+kb1))*E1 - ((dmE*T1_)/(T1_+kd))*E1 - dE1*E1 - gE*((T1+T1_)/(T1+T1_+kg))*E1 + nE*((aE*VI)/(VI+kV))*E2
		ddE2 <- gE*((T1+T1_)/(T1+T1_+kg))*E1 + ((bE2*kb2)/(E2+kb2))*E2 - dE2*E2 - ((aE*VI)/(VI+kV))*E2
		list(c(ddT1, ddT1_, ddT2, ddT2_, ddVI, ddE1, ddE2))
	})
}

library(deSolve)
X <- ode(X0, seq(0,100,by=0.01), fn_ode, params)

colnames(X) <- c('time',paste('x',1:7,sep=''))  # convert back to old nomenclature
